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ABSTRACT_ 


A numerical model which utilizes the isobaric vorticity equation 
is developed and applied to cumulus-scale data. The model, together 
with a modified version of the cumulus convection model of Weinstein 
and Davis, is applied to data obtained from the National Severe Storms 
Laboratory in Norman, Oklahoma. The calculations yield real time pres 
dictions for height, temperature and olnenee humidity at seven pres= 
sure levels, which are then used as input to the cumulus convection 
model to obtain vertical profiles of various parameters at specified 
grid points. 

Some results of the calculations are presented along with sug- 
gestions for further testing and improvement. The results indicate 
that further modifications to the approach used are necessary in order 
to provide more accurate forecasts. Values of the individual terms 
in the vorticity equation are presented as computed from the observed 


mesoscale data. 
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1. INTRODUCTION 


The National Severe Storms Laboratory (NSSL) located at Norman, 
Oklahoma, has been gathering vast amounts of data during the thunder- 
storm season for the past several years. Many studies have been 
conducted with the aid of this data (for example see Barnes, 1968, 
Newton and Fankhauser, 1964 and Hammond, 1967), however, after an 
extensive literature search, it appears that no research has been 
done, or at least reported, involving a numerical modelling approach 
whereby the distributions of various parameters were analyzed and their 
values numerically forecast for a future time. Such an approach, if 
successful, could give distinct advantages in predicting, for example, 
where hail might be expected, or where the most favorable conditions 
for tornadoes might be located. In addition, during severe convective 
activity radiosonde data often terminates near mid-tropospheric levels. 
Consequently, a model which could reasonably predict the sounding 
above termination level, when given a complete sounding from the 
previous observation time, would aid greatly in later research efforts. 

Figure 1 shows the region used for this study with a height-field 
analysis and the reported winds at the nine radiosonde stations in 
the NSSL upper-air network (see Appendix A for a description of the 
NSSL upper-air network). Assuming the data to be free of errors, 
the non-geostrophic nature of the flow on this scale is immediately 
evident, and the figure emphasizes the need for intensive research 
into the dynamic mechanisms involved under such circumstances. In 
an effort to evaluate the possibilities of numerical weather pre- 
diction for such small-scale data, this paper uses a numerical 


prediction model based on the complete isobaric vorticity equation. 
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Figure 1: Height analysis with wind directions and speecs (m sec |) 
at nine radiosonde stations of NSSL network for thy 300-mb levee. 


1700Z on 30 May 1967 which illustrates the non-geostrophic nature o- 
the data. 


14 


The resulting forecasts are then used as input to a cumulus convection 
model developed by Weinstein and DERI = (1967) to provide an aid to 
severe storm forecasting. 

The parameterized cumulus convection model by Weinstein and Davis 
(1967) calculates cloud top, profiles of vertical velocity, tempera- 
ture excess (that is, the temperature of a parcel of air in the cloud 
minus the temperature of a parcel of air in the environment), hydro- 
meteor water content and updraft radius of the cloud as well as rain- 
fall amount and duration. The principles of the above calculations 
are based on the entrainment principles of Stommel (1947) and the 
results of Kessler (1967). With a standard radiosonde sounding at 
mandatory and significant levels beginning at cloud base, the W.D. 
model transforms the sounding to one which is spaced in equal incre- 
ments of height. After the user specifies certain boundary conditions 
such as ice nucleation temperature, radius of the base of the cloud 
and conversion and collection rates (Kessler, 1967), the model 
specifies cloud extent. One use of this model is to compare the 
effect of seeding a given cloud or allowing it to develop naturally. 
This is accomplished by varying the ice nucleation temperature. 

The numerical model presented in this paper predicts heights, 
temperatures and relative humidities at individual grid points for 
seven isobaric levels (950 mb, 850 mb, 700 mb, 500 mb, 300 mb, 200 
mb, and 150 mb) so that predicted information can be used as input to 


a modified subroutine form of the W.D. model to yield the vertical 


‘The model developed by Weinstein and Davis has been referred to as 
the Pennsylvania State model. In this paper, it will be abbreviated 
as the W.D. model. 
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profiles previously discussed, under the assumption that cloud base 
is at 850 mb. However, the present version, as described in this 
paper, does not provide for vertical coupling of the seven isobaric 
levels after the first time step. 

Predictions were made at 1.5-hour intervals utilizing the IBM 
0S/360 (MVT versions 15, 16) computer. The 1.5-hour prediction 
interval was selected to correspond to the times when verification 
data were available from the NSSL radiosonde network. 

The grid used for the computations is a 24 x 24 X-Y grid oriented 
000-180, 090-270, having a mean latitude of 35N. The grid distance 
is 5 nautical miles with the grid arranged so that the nine upper- 
air stations are either located exactly on a grid point or are so 
close to one that they may be approximated to be at the nearest grid 
point. The farthest distance away from a grid point for a particular 
station was 2 nm. They were assumed to be located at the nearest 
grid point to eliminate the need for an interpolation scheme in the 
objective analysis used for initialization of the data as described 


in Section 4. 
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2. MODIFICATIONS TO THE W.D. MODEL 


While working with and studying the W.D. model, it became apparent 
that several modifications could be made to it that would provide more 
accuracy. Although most of the modifications were of a minor nature, 
three of the modifications are considered to be significant improve- 
ments over the original model. These three modifications will be 
presented in detail while the remainder of the changes will be sum- 
marized at the end of this section since they are considered to be 
self-explanatory. 


Kessler (1967) proposed the equation 


aE = = ae = k, (m-a) (Zeta) 


to model the conversion of cloud water to hydrometeor water. Here, 
: — é sll 

ky is an empirical constant which is equal to .001 sec when m>a 
and is zero for ma. The quantity "a'' is defined as a threshold 
value and refers to the cloud water content. If the cloud water 
content is greater than this threshold value, then conversion of 
cloud water to hydrometeor water will occur. If cloud water amounts 
less than the threshold value are present, then conversion of cloud 


water to hydrometeor water will not take place. This equation is 


incorporated into the W.D. model (Weinstein and Davis, 1967) in the 


form 
on = « we —y re (Q - ') (2-1b) 
dt dt Ley 

where k, is defined in a similar manner as it was in (2-la). Kessler 


1 


(1967) indicated that measurements have been performed which indicate 


=o 
that cloud water amounts greater than l emm ” are usually associated 
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with the production of precipitation. Furthermore, he arbitrarily 
suggested the value of 0.5 gm aes for the threshold value. Weinstein 
and Davis (1967) chose a' = 0.5 gm Ket for their threshold value. 

It is proposed to consider this threshold value to be a function 
of density. This approach thereby associates a different threshold 
value with each level by making it a function of height. Using this 
approach, the conversion of cloud water content to precipitation in 


the W.D. model has been modified to 


dQ 
= = k,(Q,-a/p) (2-2) 


where 9 = air density. By defining ky as it was in (2-la), conversion 
of cloud water to precipitation is modelled to occur when oo i=) Bi 
Computations have been made to determine the variation of the 
quantity a/p. These values show a maximum of approximately .00137 
gm oni in the upper levels and a minimum value of .00076 gm at 
near cloud base. 
The second significant modification is the inclusion of a linear 
variation of the latent heats of vaporization (L.) and fusion (L,) 
with temperature. Since 


Lees Leni eee 


and Le may be considered constant while L, is a function of temperature, 


i 
then L is also a function of temperature. The values of L, (Listas 
1951) are 48.6 cal aS at -50C and 779.7 cad at OC. The latent 
heat of vaporization varies from 629.3 cal aap at -506 to 5975 


cal se at OC. Since Le varies by only 1 cal die from OC to -40C, 


its value was taken as constant and equal to 678 cal a The 
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variation of only one of the latent heats on the right side of equation 
(2-3) needs to be computed since the value of the second one can be 
obtained by subtraction. 

In order to approximate the variation of the latent heat of vapor- 
ization as a linear function of temperature, it was necessary to deter- 
mine an optimum value for the difference between the specific heats of 


water and ice. [Integrating the formula 


ie ne go on 


and applying the mean value theorem for integrals yields 


L = ay + or = Cc) (T-T,). (2-5) 


Utilizing the values of L (List, 1951) for various temperatures (T) 
and with T = -40C, the value of (Co - c,) that satisfies equation 
(2-5) was computed for each temperature at 10 degree intervals from 
-40C to +40C. A plot of cay ce versus temperature is shown in 
Figure 2, Based upon the results it was decided to select the value 
of -.62340 as an "optimum" value for (Ce ace) and consider it to be 
a constant. This value is considered to be a representative value 
between -40C and OC, and corresponds exactly to the value at -10C. 
The final significant modification involves the values used for 
the constants when computing saturation vapor pressure (Weinstein 
and Davis, 1967). Integrating the Clausius-Clapeyron equation one 


obtains the following for the saturation vapor pressure over water: 


ego are aor muon cere 2 ee aes 
2 RT R 
ol suo veey 2 py enw oo 
MZ R SO’ (226) 


LN, 


( 


a 
a 8s 


M ad 
-oinjertodwaqy yqVImM ( 9 - 


O ee 


“ ad 


(9 ) eanjeirsdway 
Oo 


MeyORUOTIe Tae) 


27 a1nsiy 
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Simplifying yields 


x 
In e = - — - B’ im T + C' 
Ss ls 
where 
Pio st ee = eee) | 
A' = ae = 6958.9262 
e (ee - ee 
Bi = - ——PY = 5.65567 
R 
and 
E 506 a ora 
1a Slee n = 
C R im (Ce . (1+1n tT) | + In c= 59.01383 
———-—— -lo -1 
when Rh = Zeselok ee - C ) = -0.62340 cal gm K 
O pv W 
and Li. = 61,7 cal a ae 
vo 
Similarly, when freezing occurs one obtains 
Ss " ” 
In e. 7 + C (2-7) 
where aA br 
A" = Fe = 6151.0205 
: _ =i 
and C" = 7 + In ee 24.3277 when TS = 233.16K and L. = 678 cal gm 
O 


The above equations yield the saturation vapor pressure in millibars. 
Obviously the number of significant digits shown above exceeds the 
accuracy of some of the assumptions made to arrive at the values shown. 
Rounding to five significant figures could be done with no loss of 
accuracy, however they have been retained here as they were computed 
on the IBM O0S/360 computer. 

There were several minor changes made to the W.D. model. These 
included changing the value being used for the acceleration of gravity 
Eom 9eev ™ sea to a more exact form obtained by computing the value 
from the formula (List, 1951) expressing g as a function of latitude, 


altering the values used for the constants for conversion from degrees 


Zu 


Celsius to degrees Kelvin, conversion from meters to feet, and the 
value of the gas constant for dry air. In addition, the section of 

the program dealing with corrections due to the wind shear were com- 
pletely eliminated which consequently allowed a subroutine (subroutine 
CHEQ) to be eliminated. Other minor changes included removing factors 
of 10 so that all pressure values are carried in units of millibars, 
introducing symbols wherever possible to simplify future modification 
and testing, changing the estimate of the updraft velocity at cloud 
base from 2.0 to 0.5 m ee and eduneraiae certain format statements 
to improve the headings on the output. Finally, many soundings were 
run through the final modified version of the W.D. model to ensure that 
all changes had been properly programmed and to ensure that the results 


from the model were reasonable. 


Ze 


3. MATHEMATICAL DEVELOPMENT OF A FORECAST MODEL 


This section describes the model used for this research. Although 
it is not consistent in many respects, it did provide a means by which 
forecasts of height at 1.5-hour eee oul be obtained for input 
to the W.D. model as well as provide a method by which the values of 
the various terms in the prediction equation could be ascertained. A 
description at the end of this section describes an attempt to set up 
the prediction equation in a consistent manner after the original 
model was running properly, however, the more consistent method failed 


to provide reasonable forecasts. 


A. THE VORTICITY EQUATION 
Following the procedures given by Haltiner and Martin (1957), the 


isobaric vorticity equation may be written as 


: yy 4 (2 BBW 
CO) ee Ga ON Gal ce et oe Se Gaia) 


a la 
ck 


Since the entire grid used for this research covers only 2 degrees 
of latitude, the Coriolis parameter can be considered constant, 
taking on the value it would have at the mean latitude of the grid. 
Using this approximation and expanding the left side of (3-la) yields 


eS Bae we ov. woe 
Nea (w$ eV EVE SE (Ct+f Vv Wot Ses aes a: (3-1b) 


In order to use the W.D. model at 1.5-hour intervals, height of the 
base of the clouds is needed from which computations can be started. 
To obtain these heights, (3-lb) will be used as a height forecasting 
equation analogous to methods applied to large scale motions in the 
atmosphere by applying the balance equation. An alternative method 
without using the balance equation for initialization was tried and 


1s discussed at the end of this section. 
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Since the application of (3-1b) to height forecasting will be done 
to obtain information regarding the orders of magnitude of the various 
terms and to obtain heights which can be used as cloud base information 
for input to the W.D. model, no attempt was made to set up a consistent 
set of equations. To obtain an equation which could be used for height 
forecasting and evaluation of the various terms with mesoscale data, 
the relative vorticity on the left side of (3-1b) and in the diver- 
gence term was replaced by Vy. In addition, only the rotational wind 


components were used for the horizontal advection term, that is, 


u = - 3 and v - 2 (3-2) 


This was done because the present version has no provision for pre- 
dicting the divergent component of the wind CW). Future modifications 
should provide the means by which we can be used in addition to fore- 
casts of the rotational component CW). The vorticity in the vertical 
advection term as well as divergence and the u and v components of the 
wind in the twisting term were computed based on the observed winds 

at the initial time and held constant for 1.5 hours. Keeping values 
constant in the - terms greatly Staplliweties programming considerations 
since 3-minute time steps were used and to update these values at each 
time step would involve working with each level for only one time step 
before proceding to the next level. Unfortunately, holding C in the 
vertical advection term and the u and v components of the wind in the 
twisting term constant for 1.5 hours causes the seven levels to lack 
vertical coupling. Future modifications should attempt to correct this 
deficiency. After making these substitutions and assumptions, the 


resulting prediction equation becomes 


Sy) = - (Se - HK WK, ye BW _ BA BG. 
Sa y) = - (os So + ot ll Voy + fv: WW + Ro ayo (323) 
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Assuming continuity in t, the term on the left side of (3-3) may 
| 2 dy | | 

be written as V sya: and the equation then becomes a Poisson-type 
equation which can be solved by relaxation procedures as described by 
Haltiner and Martin (1957). The forecasting procedure, being one of 
successive iterations, will thus allow all terms involving ¥ in (3-3) 
to vary at each time step. 

To obtain a stream field for (3-3), consider the balance equation 


(Gharney, 1955): 


Vy += (WHVE) - 27 (ee Ly. ais nud (3-4) 


Since the Coriolis parameter is to be considered constant over the 
grid, the second term of (3-4) drops out. 
As discussed by Charney (1955) and Bolin (1955), equation (3-4) 


will maintain its elliptic character only if 


INO 
Da 


Vi4+5> 0, (Ee) 
that is, only if the geostrophic relative vorticity is greater than 
-£/2. As shown in Appendix B, heights would have to be read at least 
to the nearest centimeter on this size grid in order to satisfy the 
requirements of (3-5) and still allow negative geostrophic relative 


vorticities. Consequently, (3-4) was forced to be elliptic by 


writing the non-linear terms in the form 


ou ov 
EL). ~ Ox Oy ? encey 
and 
oe Fy _ Ov ou (3-6b) 
See ae ~ Ox dy 


Making this substitution into (3-4) and taking into account that 
the Coriolis parameter is to be considered a constant, (3-4) becomes 


vy = By, - = ae Gt vis ap (3-7) 
m 


m 


ZS 


where J represents the Jacobian operator, and u and v are the components 
of the observed winds. Equation (3-7) is then an elliptic equation as 
it stands and does not require that the geostrophic relative vorticity 
be greater than -f/2. The field can be initialized by the approxi- 


mation = ee, and then equation (3-7) can be relaxed to convergence 


m 
before using (3-3) to predict the ¥ fields at a later time. Similarly, 


once the ¥ fields have been forecast ahead, the height fields can be 
recovered by making the approximation z = = and then utilizing (3-7) 
imeéhe relaxation preeeduge untulconvergence asmateaimed: 

Convergence was defined as occurring when (3-7) balanced within 
O.5 meters. When converting from height fields to stream fields 
this requires an epsilon of approximately 58612 ma cca whereas 
epsilon becomes 0.5 m when converting from stream fields to height 
fields. These values arise from the approximation ¥ = (g/f )z. In 
retrospect, the above convergence criteria may be too large for this 
scale of motion. In the future, more stringent criteria (such as 0.1 
m) should be invoked. 

As mentioned at the beginning of the section, attempts were made to 
correct the deficiencies of the prediction equation. This involved 
utilizing the exact relation v4 = C as the method of obtaining the 
initial J field instead of using the balance equation and to substitute 
v4 for ¢ in the horizontal advection term. In addition, the conver- 
gence criteria was changed to 10000 AC pcm in units of }, when con- 
verting from € to ¥. 

The balance equation was used to recover height fields from the 
forecast fields of ¥ with convergence criteria of 0.1m. Two methods 


of initializing the | field were tried; they were = ra and ¥ =e 


m m 
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where Zz represents the standard atmosphere value of the height of the 
respective isobaric surfaces. When using z for the initialization, 
computations were carried out without difficulty, however the result- 
ing forecast height fields were very flat as might have been anticipated. 
However, the forecast fields were in error by as much as 300 m at some 
levels. When using z for initialization, instability resulted in the 
horizontal advection term which caused an exponential growth of this 
term to values too large for the computer to handle after only 12-15 
time steps (36-45 minutes ahead in time). The v4 part of the hori- 
zontal advection term is believed to be the cause of the instability 
but this was not investigated in detail. However, one attempt to 
investigate the instability was made by shortening the time step to 

one minute. This delayed the instability somewhat which would indicate 


the problems were caused by computational instability. 


B. VERTICAL MOTION COMPUTATIONS 

In order to solve (3-3), vertical motion values are required for 
the vertical advection of relative vorticity and for the twisting term. 
Vertical motions for each of the seven levels were computed by means 


of the continuity equation in the form 


ee 
paras (3-8) 


where the bar symbol denotes a vertical average through the layer. 
This method of computation utilizes the measured wind fields and does 
not allow the w values to be updated for forecast times since the 
present model has no provisions for obtaining the divergent component 
of the wind. The vertical boundary conditions imposed were 

©1000 = one) (3-9a) 
and Ww = 0.0 (3-9b) 


1ESN8 


reas YAY 


which also imply that V° W = 
1000 p = 950 


Vertical motion computations were made by working downward from 
150 mb through 300 mb and upward from 1000 mb through 700 mb. Vertical 
motion at 500 mb was then computed as the average of the values at 300 
mb and 700 mb. It is recognized that using a backward or forward finite 
difference form for could at times lead to computational difficulties, 
however the method used here does allow the vertical motion at a given 
level to be affected by not only the divergence at that level and the 
level above, but also takes into consideration the divergence field at 
the level below. 

Computations for the various levels were as follows (symbols such 
as W150 and W950 were chosen to conform to the symbols used in the 
computer program (subroutine VERTMO) shown in Appendix C): 

(1) 150 mb 

W150=0.0 
(2) 200 mb (Ap>0) 


W175=W150-ApV> Wo o9_ 159 


W225=W150-ApV* Wa59_ 150 


W200=(W175+W225) /2 
(3) 300 mb (Ap>0) 


W2 50=W200-APV> Wy 50 


W350=W200-ApV> We 909 
W300=(W2 50+W350)/2 
(4) 1000 mb 


wicoc=c.9 
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(5) 950 mb (Ap<0) 


———— a 


W925=W1000-APV> Ween sq 


W975=W1000-APV> Woe ig 
W9 50=(W925+W975)/2 

(6) 850 mb (Ap<0) 
W825=W950-ApV* Wo 04 _ 5, 

W9I00=W950-APV> Ween gc, 
W850=(2 (W825)+W900)/3 

(7) 700 mb (Ap<0) 


W675=W850-ApV. Wenn gen 


W725=W850-ApV- V5 00-950 
W700=(W675+W725)/2 

(8) 500 mb 
W500=(W300+W700)/2 


The values obtained for vertical motion by this method are discussed 


iimseceLOnN 5. 


C. TEMPERATURE AND RELATIVE HUMIDITY FORECAST SCHEME 

In order to utilize the W.D. model as an aid in severe storm fore- 
casting, forecasted values of temperature and relative humidity are 
necessary in sounding form as input to the model. 

Although a model designed for operational use should definitely 
be internally consistent, such as the vertical consistency of tempera- 
ture fields with the geopotential fields, this study involved deter- 
mining some of the problems that would develop when applying numerical 
forecasting techniques to mesoscale data. Consequently, the scheme 
for obtaining temperature forecasts which follows is not dynamically 


consistent with what has preceded, but it does provide a method of 
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obtaining a temperature-field forecast which could be used as input 
to the W.D. model. As discussed in Section 6, future studies should 
improve on this method. 

Four different methods of forecasting temperature were tested and 
the results were then compared to verifying data. The method which 
verified best was selected as the forecast method for temperature, 


which was simply 


OT _ | 
<a Ly) 2a Eq 0) 


where the values for the wind components were taken as the measured 
winds at observation time and were not altered during the forecast 
interval. 

Similar testing of four prediction methods was also done for 


relative humidity. The results led to 


OR, 
ae yy \WeVR, (3-11) 


as the prediction equation for relative humidity where the wind 
components were again taken as the measured winds at observation 
time and were held constant for the forecast interval. 
Temperature and relative humidity prediction calculations by 
equations (3-10) and (3-11) were performed in subroutine PROG, a 
listing of which may be found in Appendix C. Results of the com- 


putations, when compared to verifying data, are found in Section 5. 


D. BOUNDARY CONDITIONS 
To initialize the | field, the approximation } ~ gz/f_ was used 
as discussed in Section 3A. When the relaxation was performed, the 


outside row or column of the grid could not be modified due to finite 
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differencing. As a result, the boundary values of ¥ were held at the 
value obtained by the initialization. The prognostic equation (3-3) 
was utilized from x = 3, 22 and y = 3, 22. After initializing the 
on field to zero, it was then relaxed for the first time step. This 
relaxation was performed from x = 3, 22 and y = 3, 22 as mentioned 
above and the outer two rows or columns were maintained at their 
original values of zero. By this means, actual forecast values of 

$ were obtained only from x = 3, 22 and y = 3, 22. After the first 
time step, the previous relaxed values a = were used as the initial 
estimate for the new time step. 

After obtaining the 1.5-hour forecast for the ¥ field, the con- 
version back to z fields by means of (3-7) was accomplished by making 
the initial approximation z ~ ff t/g only for the points where the ff 
field had been forecast. This then leaves the outer two rows or 
columns without forecasted height values. Before relaxing the z 
field, however, the values at these boundary points were set by 
performing a linear extrapolation from the fourth row or column 


OWEward as follows: 


4 
3 
2 Zo = 2, - (2), - Z) 
1 Z, = Zy - (z., - Z5) 


Once this was accomplished at all of the grid points in the two 
boundary rows or columns, the z field was relaxed by means of the 
batliance equation from x = 2, 23 and y = 2, 23, 

Horizontal boundary values were not a problem with the vertical 


motion fields since they were computed from x = 2, 23 and y = 2, 23. 


on 


As mentioned above, the prognostic equation was utilized from x = 3, 22 
and y = 3, 22 so that there were computed values of omega at x = 2, 
x = 23, y = 2 and y = 23 to use for the terns © and 2. As stated 
in (3-9), vertical boundary conditions for w were that the vertical 
motion would vanish at 150 mb and 1000 mb, which then implies that 
the divergence at 1000 mb is equal to the negative of the sum of the 
divergences from 950 mb to 150 mb. 

As will be discussed in Section 4, the initial data fields of 
height, temperature, relative humidity, wind direction and wind 
speed were analyzed on the computer by an objective scheme. The 
analysis method included finite differencing methods such that the 
outer row or column was not analyzed. These boundary values were 
set to the same value as the next interior row or column. 

The same linear extrapolation technique that was used for the 
height fields was also used for the Bouncer, conditions of the fore- 
casted values of temperature and relative humidity, but it was only 


necessary to extrapolate from the second row or column out to the 
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E. TIME STEPS 

In the forecasting of the new } fields by (3-3), a 3-minute time 
step was used, thus requiring 30 time steps to obtain a 1.5-hour fore- 
cast. The method used was the leap-frog scheme described by Haltiner 
(upbl.), with a forward time step used for the first step. To 
eliminate the possibility of obtaining two separate solutions by thie 


leap-frog method, the two solutions were averaged every 10 time steps 


SZ 


(30 minutes) and then one pass through the averaged field was made 
utilizing a Laplacian smoother with a .04 coefficient for smoothing. 
The leap-frog forecast ne thod was then continued. 

The forecasting of temperature and relative humidity by the simple 
advection schemes described in Section 3C does not involve any 
iterative type computations since the wind components were held 
constant for the forecast interval. Consequently, a single 90-minute 
forward time step was used. However, future modifications, which 
should include provisions for obtaining rotational and divergent com- 
ponents of the wind, should also update the wind field with each time 
step so that forecast values of the wind field could be used for 


temperature and relative humidity advection. 
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4. INITIALIZATION OF THE MODEL 


In order to commence computations with the previously discussed 
model, 35 fields of data are needed (7 levels each of height, tempera- 
ture, relative humidity, wind direction and wind speed). Such a 
volume of data would be difficult to obtain by manual means for 
real-time operational use of a model. Therefore, a more realistic 
approach is to use an objective analysis scheme by the computer. 

Although many complicated and very sophisticated analysis schemes 
can be and have been developed (for example see Hughes, 1967), this 
was not the purpose of the research reported in this paper. Con- 
sequently, a simple objective analysis scheme was devised which uses 
the Laplacian operator and an eight-point averaging technique in 
alternating steps. Comparisons of computer analyses and hand analyses 
showed the patterns to be very similar. The major difference between 
the two analysis methods was the gradient in the vicinity of closed 
centers. The hand analyses tended to subjectively spread the gradient 
out over a much larger area whereas this was not so with the computer 
analyses. Which analysis is more nearly correct is merely a matter 
of subjective judgement. With only nine data points within a 14,400 
square mile area during conditions of severe storms, this question 
will probably not be answered until the data network becomes more 
detailed. Nevertheless, the computer Aaa brace. were judged to yield 
adequate data fields for testing of the computer program. 

In order to start the analysis a first guess was needed at each 
grid point. This was accomplished by dividing the grid into nine 


"regions of influence" (see Figure 3), each region containing one 
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of the nine radiosonde stations. The value of a given parameter at 
each grid point within a region was then set equal to the known value 
at the radiosonde station. 

With the initial guess completed, the analysis was started. The 
first step was to alter the value at the grid points (except for the 
nine data points) by setting the value at that point equal to the 


average of the surrounding eight points according to the following 


diagram: 
4 3 2 
5 0 ly 


gh 2 ere, Sy Cl 


0 8 
After making one pass through the grid in this manner, the Laplacian- 


type scheme was used such that 


P Pa +o (P,P Pa tee ee (G82) 


0 0 3 5 i, 


where the value of 6 O.l1 was arbitrarily selected. After making 
ten passes through the grid in this manner, the eight-point averaging 
technique was again employed. The entire procedure was repeated in 
such a manner that the eight-point average was employed 11 times while 
the Laplacian scheme was utilized 100 times. The resulting analysis, 
after boundary conditions were set, was ene analysis used in the pre- 
diction model. The boundary condition imposed was that the value at 

a grid point on an outside row or column was equal to the value at the 


adjacent grid point on the next row or column toward the center of the 


Cunt ¢ 
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Figures 4a and 5a show hand analyses for two different levels and 
Figures 4b and 5b show the corresponding computer objective analyses. 
A copy of the program for the objective analysis scheme is included in 


Appendix C, 
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Figure 4a: Hand analysis for 950 mb at 1700Z on 30 May 1967. 
Contours every 5 meters. 
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Figure 4b: Computer ol-jective analysis for 951 
sO May 1967... “Contours Jee Cer Ss. 
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Figure 5a: Hand analysis for 500 mb at 17002 on 30 May 1967. 
Contours every 10 meters. 
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Figure 5b: Computer objective analysis for 500 mb at 1/700Z on 
30 May 1967. Contours every 10 meters. 


41 


5. RESULTS OF THE COMPUTATIONS 


Two sets of data have been utilized with the previously described 
model. The times of the two sets of data are 1/00Z, 30 May 1967 and 
2130Z, 30 May 19G7: 

In synoptic-scale models, the terms in equation (3-la) involving 
vertical motion and divergence have often been neglected since they 
may be considered relatively unimportant by scale analysis (Haltiner, 
upbl.) in the large scale motions of the atmosphere. However, when 
dealing with cumulus-scale data, this may not be true. In order to 
evaluate the relative importance of the various terms in the model 
throughout the entire grid, computations of the orders of magnitude 
of the various individual terms of equation (3-3) have been performed 
for both the 2130Z and 1700Z sets of observed data for the seven levels 
in the vertical described in Section 1. Tables 2 and 3 show the mini- 
mum and maximum values of each term based on the observed data, not 
on the forecast values. Although these tables only show the minimum 
and maximum values, they do point out the fact that the divergence 
term is very important. There were many instances when all five 
terms were of the same order of magnitude, that is, when each term 
was just as important as each of the other terms. Table 4 shows a 
sample of the values of each of the five terms as extracted from the 
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Out of 4800 values tested, Se Gp and se Oe were of the same sign 
2584 times, or about 54% of the time. Therefore, in the cases tested, 
they were additive more often than not which is further justification 


that they should be retained for small-scale data. 
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Vertical motion computations yielded considerably different values 
for the 1700Z data and the 2130Z data. Table 5 shows the maximum up- 
ward and downward vertical motions for each isobaric level at each of 
the times. Although the values in Table 5 may, at first, seem somewhat 
large, a transformation by the formula 

aw -p Bw (5-3) 
where 0 is replaced by p/RT from the equation of state to yield 
we - £ ow, (S=) 
RT 
shows that a value of w = 200 mb * i upward corresponds to a value 
of w .86 m sae upward (at 500 mb) which is certainly a reasonable 
value for cumulus convection. Although values this large were only 
present at a few points, further research is needed to determine what 
reasonable values are for this scale of motion. In general, the 
vertical motion fields showed fairly strong upward vertical motions 
throughout the grid at 1/700Z and correlated fairly well with the 
squall-line that was present at that time. There were a few areas 
in which the correlation was not too good but this may possibly be 
due to the fact that the nine radiosondes were not all launched 
simultaneously and a difference of 15 minutes between launch times 
could possibly result in the bad correlation regions. Although the 
radar pictures indicated that dissipation of the squall-line had 
begun at 2130Z by the fact that isolated echoes were present rather 
than a large mass of clouds, the vertical motion fields did not show 
extensive regions of downward motions except at 300 mb. The 21302 
values were, however, considerably less in intensity than they were 


at 1700Z which might have been anticipated from the results shown in 


Tab le wae 
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The significance of the much larger values of vertical motion at 
1700Z as compared to 2130Z was not investigated thoroughly, however 
larger values of divergence were evident (as would be expected) at 
1700Z. This was particularly true at 200 mb, and this would tend to 
be propagated downward to 300 mb and 500 mb by the computation procedure. 

The simple method of forecasting temperature resulted in values 
which verified slightly better than a persistence forecast. The 
average absolute error (not considering whether the error was an 
overestimate or an underestimate) for 97 persistence temperature 
forecasts was 1.8C while for 105 50% temperature advection forecasts 
the average absolute error was 1.4C. The reason for the different 
number of verifications is due to the fact that some of the 1/00Z and 
2130Z radiosonde data terminated before reaching 150 mb. Verifications 
were performed only at the nine radiosonde stations, and only at each 
of the seven discrete isobaric levels. 

The simple method of forecasting relative humidity resulted in 
values which verified essentially the same as a persistence forecast. 
The average absolute error for 9/7 persistence relative humidity fore- 
casts was 18.0% while for 105 50% relative humidity advection forecasts 
the average absolute error was 18.2%. 

The results of the height forecasts did not verify as well as was 
expected, This is undoubtedly due in part to the fact that additional 
refinements are necessary in the prediction method to obtain a con- 
sistent set of equations. Table 6 compares the average absolute 
errors of persistence and model 1.5-hour forecasts for both the 1/00Z 
and the 2130Z sets of data. As can be seen in the table, the model 


did perform slightly better than persistence 50% of the time. 
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As previously mentioned, the W.D. model yields profiles of vertical 
velocity, hydrometeor water content and temperature excess. [In addition, 
total rain expected from the cloud, duration of the precipitation, up- 
draft area, and a profile of updraft radius are included in the results 
of the computations. Figure 6 shows a sample of the graphical output 
of the W.D. model. These results are printed out in graphical form 
by the two subroutines (GRPHCL and GRAPH) which were included with the 
W.D. model (Weinstein and Davis, 1967). The horizontal axis (shown in 
Figure 6 at 1/2 scale) shows the scale to be used for the various para- 
meters and the vertical axis on the left shows the corresponding 
heights in meters, with the bottom figure (in this case 1642) being 
the first increment (200 meters) above cloud base. In the actual 
computer output of Figure 6, the lines are represented by letters, 
such that the vertical velocity curve is a series of the letter "Ww", 
hydrometeor water content of the letter "Q'"', temperature excess of 
the letter "T", updraft radius of the letter "R'', and the zero axis 
for the temperature excess curve is represented by the symbol "$", 

The W.D. model has a provision for setting a constant value for 
the updraft radius or for allowing it to vary within a cloud. How- 
ever, for this research the updraft radius was always considered to 
be constant and equal to 1 km. The profile shown in Figure 6 was 
computed based upon the 1830Z,30 May 196/ forecasted values of 
temperature, relative humidity and height of the cloud base, and 
by assuming that cloud base was at 850 mb for station CHK. 

Many soundings were used as input to the W.D. model, and in many 
cases the values of the temperature excess determined by the model 


seemed to be too large. Values as high as 4 or 5 degrees Celsius 
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occurred quite often, with an occasional value of 7 or 8 degrees 
Celsius. This intuitively seems to be too high but additional 
observations are needed to determine if such values actually exist 
within cumulus clouds. In addition, apparent discontinuities in 
temperature excess, as evidenced near 8000 meters on Figure 6, are 

quite common in the graphical output of the model. These large tem- 
perature excess values and temperature excess discontinuities occurred 
not only with forecast soundings but also with actual observed soundings 
from the NSSL data. 

Figure 6 is based upon an ice nucleation temperature of -25C which 
is assumed to be representative of natural conditions in the atmos- 
phere. The W.D. model allows for changing the ice nucleation 
temperature to a value of, say, -6C to represent the conditions that 
might occur if cloud seeding were to be undertaken. However, this 
is presently not included in the subroutine version of the W.D. 
model (subroutine CLOUDS) since the reason for using it as part of 


a prediction model was to aid in severe storm forecasting. 
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6. SUGGESTIONS FOR FUTURE DEVELOPMENTS AND IMPROVEMENTS 


The first modification to the previously discussed model, which 
became evident during testing of the model, would be to expand the 
horizontal grid dimensions to at least 32 x 32. Figure 3 illustrates 
that four of the nine grid points lie within the first two rows of 
the boundaries. This more or less wesateet application of the prognostic 
equations to these grid points due to finite differencing. It is esti- 
mated that this improvement would increase the storage requirements of 
the program by about 50%, and since the program is already quite large, 
modifications would have to be made to allow for the sharing of some 
storage locations wherever possible. 

At present the model lacks the proper provisions for predicting 
the complete wind field and for updating it after each time step. 

A consistent scheme for updating the wind field should be included 
such that the predicted winds could be recovered every time step. 
In addition, this would allow predictions further ahead in time as 
well as allow for updating the vertical motion fields with new 
divergence fields. 

As discussed in Section 3A, the present version of the forecast 
model does not provide for vertical coupling of the seven isobaric 
levels since the vorticity in the vertical advection term and the u 
and v components of the wind in the twisting term were held constant 
for the 1.5-hour forecast interval. To correct this deficiency would 
involve performing one time step at a given level before proceding to 
the next level. Also, additional storage would be required to save the 
field which Operates on for one time step back so that vertical 


derivatives could be taken at each level on fields which all correspond 
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to the same time. As far as programming considerations are concerned, 
this essentially involves nesting the vertical level DO LOOP inside 
of the time step DO LOOP. 

When the model has been modified to allow for vertical coupling of 
the seven isobaric levels, additional levels in the vertical could be 
included so that vertical layer depths would not exceed 100 mb. This 
would require perhaps 10 or Il levels which would increase the storage 
requirements of the computer considerably, however, this would allow 
a much finer detail for the forecasted soundings. [In addition, it 
should aid in the delineation of the vertical velocities at each 
level. 

A consistent dynamical method should be included for the prediction 
of temperature and relative humidity. The present method is more 
empirical than dynamical and better verifications should be obtained 
if a more realistic approach is devised. As an alternative, a better 
COs Gali hee 50%) might be tested, particularly in the case of 
relative humidity where the errors were negative much more often 
than they were positive. 

In subroutine SETUP (see Appendix C), cloud base is assumed to be 
at 850 mb for entry into subroutine CLOUDS (see Appendix C). A routine 
could be developed to be included in SETUP such that the lifting con- 
densation level, mixing condensation level or convective condensation 
level would be computed as the cloud base before entering the W.D. 
model (subroutine CLOUDS). Comparison of results obtained for cloud 
base at 950 mb and 850 mb indicate that the cloud base and the associ- 
ated temperature and relative humidity are very critical to the results 


obtained from the W.D. model. 
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Once a complete internally consistent model is set up, calculations 
out past 1.5 hours could be tried to determine at what point such a 
model loses its skill. A few runs were made out to 4.5 hours and the 
results appeared to indicate values of height which were too large, and 
the temperature fields appeared to become unstable. 

Kessler (1967) proposed the formula for the fall speed of pre- 
Cipitation to be of the form 

V = 550939 Meee (6-1) 

The corresponding formula in the W.D. model is 


mz 


Hees) (G2) 


The units of M are gm ey whereas the units of Q. are gm an in 
order to convert M es the units of Q >» it is necessary to divide by 

the density of air, o. Weinstein and Davis (1967) have apparently 
selected a representative value for 9. It is recommended that future 
modifications of the W.D. model include the density in the computation. 
Since M = PQ» and 9 = p/RT from the equation of state, it would be a 
Simple change in the programming to use the following formula for the 
hydrometeor terminal velocity: 


PQ. wh) 
V = (580967 =) : (6-3) 


Since the duration calculation is based on the hydrometeor terminal 
velocity at cloud top, (6-3) will permit values of V which not only 
depend on the hydrometeor liquid water content but also will depend 
on the pressure and temperature at the cloud top. 

Finally, it has been suggested to the author by Dr. J. D. Mahlman 
(personal communication) that the only solution to overcome the many 


difficulties with this scale of motion would be a direct integration 


aD 


of the primitive equations. It is believed that no such research has 
been attempted for mesoscale data and this approach definitely warrants 
consideration for future attempts at numerical modelling of such scales 


of motion. 
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7. SUMMARY 


A numerical weather prediction model has been described for use 
with cumulus-scale data. Together with the W.D. model, it provides 
forecasts of height, temperature, and relative humidity at seven 
isobaric levels as well as vertical profiles of vertical velocity, 
temperature excess and hydrometeor liquid water content at 1.5-hour 
intervals. The model was tested on two sets of data and yielded fore- 
casts which were, in general, statistically about the same as a per- 
sistence forecast. This can probably be attributed in part to the 
fact that the present version has no provisions for coupling the 
seven isobaric levels nor for updating the vertical motion fields 
after each time step. Another factor is that the temperature and 
relative humidity forecast schemes were more empirical than dynamical, 
and were not consistent with the vorticity equation. This may account 
for the lack of skill in forecasting these parameters. Suggestions 
were made to correct some of these inconsistencies in the event that 
further research is conducted with the model. However, the possibility 
of applying numerical forecasting techniques to small-scale data does 
appear to have possibi ities but extensive research is needed in this 
area before useful forecasts will be obtained. 

The model requires approximately 5.5 minutes on the IBM O0S/360 
(MVT Pe eiete Le. 16) computer to yield a 1.5 hour forecast. An 
additional 14 seconds per field is required for the objective analysis 
scheme described in Section 4, thus requiring about 13.5 minutes for 
the analysis of 35 fields of data and the prediction of the afore- 


mentioned parameters for 1.5 hours. 


a7 


REFERENCES 


Barnes, S. L., 1968: On the source of thunderstorm rotation. Tech. 
Mem. ERLTIM-NSSL 38, 30 pp. 


Bolin, B., 1955: Numerical forecasting with the barotropic model. 
Tellus, 7, 27-49. 


Charney, J., 1955: The use of the primitive equations of motion in 
numerical weather prediction. Tellus, 7, 22=26. 


Haltiner, G. J., Numerical Weather Prediction, unpublished notes, 
The Naval Postgraduate School, Monterey. 


Haltiner, G.J., and F. L.. Martin, 1957. Dynamical lecind ary sme ae 
Meteorology. McGraw-Hill Book Co., Inc., New York. 


Hammond, G. R., 1967: Study of a left moving thunderstorm. Tech. 
Mem. IERTM-NSSL 31, 75 pp. 


Hughes, R. E., 1967: U. S. Naval Weather Service Computer Products 
Manual, NAVAIR 50-1G-522. 


Kessler, E., 1967: On the continuity of water substance. Tech. Mem. 
LERTM-NSSL 33, 125 pp. 


List, .R.J., Hd... 19o1: . Smithsonian Meteguokeeiicaiap less 


Smithsonian Institution. 


Newton, C. W., and J. C. Fankhauser, 1964: Movement and development 
patterns of convective storms, and forecasting the probability of 


storm passage at a given location. NSSP Report No. 22, 53 pp. 


Stommel, H-, 1947: Entrainment or air imto’a cumulus’c loud. “Joes 
Met., 4, 91-94. 


Weinstein, A. I., and L. G. Davis, 1967: A parameterized numerical 


model of cumulus convection. Report No. 11 to the NSF, (NSF 
GA-777), Dept. of Met., The Penn. State Univ., 66 pp. 


58 


APPENDIX A 


THE NSSL RADIOSONDE NETWORK 

There were 10 radiosonde stations within the NSSL region considered 
for this study which are listed below. The grid-point numbers associ- 
ated with each station are shown on Figure 3, where the stations are 


shown by their call letters. 


Location 
Call letters Lat. Long. 


Sheppard AFB, Texas 


Ringling, Oklahoma 
Altus, Oklahoma 

Fort Sill, Oklahoma 
Pauls Valley, Oklahoma 
Cordell, Oklahoma 
Chickasha, Oklahoma 
Oklahoma City, Oklahoma 
Tinker AFB, Oklahoma 


Watonga, Oklahoma 





All of these stations, except OKC, were used in this study. The 
reason for not including data from OKC is that their data corresponds 


to synoptic data times rather than at 1.5-hour intervals. 
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APPENDIX B 


ON THE ELLIPTICITY REQUIREMENT OF THE BALANCE EQUATION 
The ellipticity requirement of the balance equation (Charney, 1955 


and Bolin, 1955) was given by equation (3-5) which is 
v8 


so 


No[ Fh 


> 0. (3-5) 


Bolin (1955) stated that the data at grid points which did not satisfy 
(3-5) was modified so that (3-5) was satisfied. This section investi- 
gates the extent to which grid distance influences the requirement of 

(3-5). 


Expanding (3-5) results in 


Z 
ney eee o (Bam) 
2 Ue 
df 
m 
where 
vw = finite difference form of the Laplacian operator, 
ine = Map sdeceOn.—. hhoe 
g = gravity ~ 9.80, 
d = grid distance, 
zZ = height, 
and f = mean Coriolis parameter = mesa 


Upon substitution of the above values for each factor, (B-1l) reduces to 
2 Ze 
Vz oe .0003/7d = 0 (B-2) 
where z is in meters and d is in nautical miles. Substituting a 
value of d = 5 nm yields 
V2 > -.02175. (B-3) 


By measuring heights to the nearest meter, the only way to satisfy 


(B-3) is to have vz > 0. To satisfy (B-3) on a grid with d = 5 nm 
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requires z to be measured at least to the nearest .0l1 m, and this would 
still not guarantee meeting the requirement. It is easily shown with 
equation (B-2) that the smallest grid distance on which (3-5) can be 
satisfied, without eliminating the possibility of v2 =< -Qhand sta il 
measure z to the nearest meter, is 35 nm. 

The result from above is what prompted the substitution of the 
measured wind field into the balance equation in order to force it 
to be elliptic. Tests were made to determine with what frequency the 
balance equation ellipticity requirement was satisfied with this scale 
of data. The results show that out of 6/776 points, only 62.4% satis- 
fied the requirement. At 1/00Z the requirement was satisfied 59.4% 


of the time while at 2130Z the figure was 65.5%. 


61 


APPENDIX C 


This appendix includes a listing of the computer program used for 


the research described by the preceding sections. In order to aid 


any possible research in the future with this program, each subroutine 


will be listed with a brief explanation of the calling arguments. 


lie 


Subroutine VERTMO 


Zonal component of the horizontal wind in m sec 


U 


Meridional component of the horizontal wind in m sec 


V 


OMEGA = w = Vertical component of the wind in (x, y, p, t) 


° e -] 
coordinate system in mb sec .. 


NI = Number of grid points along x-axis. 
NJ = Number of grid points along y-axis. 
NK = Number of levels in the vertical. 


DIST = Distance between grid points in meters. 

DSI = The negative of the divergence field for 1000 mb = the 
sum of the divergence from 950 mb to 150 mb. 

Subroutine FORFUN 

VORT = Vertical component of relative vorticity. 


E = Field tor toreine function, 


FMAP = Map factor (mean). 
FBAR = Coriolis parameter (mean). 
DIST = Distance between grid points in meters. 


DT = Forward time step in seconds. 


Zonal component of the horizontal wind in m sec 


US PD 


VS PD Meridional component of the horizontal wind inm sec 
K = The level the main part of the program is at when this 


Subroutine is called. 
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IZ = Counter to determine at what time step (1-30) the main 
program is. 
; 2 ~1 
PSI = The | field inm sec . 
OMEGA = w = Vertical component of the wind in (x, y, p, t) 


coordinate system in mb sec 


Subroutine RELAXI 


M = Number from 2, NI = a number one greater than what is 

needed in DO LOOP. 

N = Same as M except 2, NJ. 

F = At field for relaxation of balance equation. 

A = field in fe cCcuee 

B = Z field in meters. 

L = Counter to keep track of number of points which have 
converged. 

EPS = Epsilon = convergence criteria. 

ALPHA = Over-relaxation constant. 


FMAP = Map factor (mean). 


Coriolis parameter (mean). 


FBAR 
DIST = Distance between grid points in meters. 
IL = Lower point for DO LOOP in x direction. 
JL = Lower point for DO LOOP in y direction. 
IM = Iteration number from main program. 


Zonal component of the horizontal wind in m sec 


Se 
tl 


<< 
il 


Meridional component of the horizontal wind in m sec. 


K = Level in the vertical from main program. 


) 
It 


Gravity. 
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Subroutine RELAK 


Same as in RELAXI except: 


EF = Forcing function: 
R = Residual field. 
ICOUNT = L. 

Il = IL. 

Jt = 22. 


Subroutine SMOOTH 
D = Field to be smoothed. 
RF = Smoothing factor (<1.0). 


IPASS = Number of smoothing passes to make through the grid. 


Jie ook Lower point for DO LOOP, 


IH, JH = Upper point for DO LOOP. 
Subroutine PROG 

1EMP)= Temperature. 

RH = Relative humidity. 


DIR = Wind directions. 


3 PD 


Wind speeds. 

USPD = Zonal component of SPD in m sec. 

VSPD = Meridional component of SPD in m sec. 
DIST = Distance between grid points in meters. 
NK = Number of levels in the vertical. 


NI = Number of grid points on x-axis. 


NJ = Number of grid points on y-axis. 


NI1 = NI-1. 
NJ1L = NJ-1. | 
OMEGA = Ww. 


Z = Height field. 
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7. Subroutine SETUP 
TEMP = Temperature. 
RH = Relative humidity. 


Z = Height field. 


IL = A number from l to 
TH = A number from 1 to 
JL = A number from l to 
JH = A number from l to 


8. Subroutine CLOUDS. 


NI, 
NI, 
NJ, 


NJ, 


usually 1. 
usually NI. 
usually 1. 


usually NJ. 


This is the modified version of the W.D. model (Weinstein and 


Davis, 1967). 


PDATA = Pressure levels. 


ZDATA 


I 


Height field. 


RHDATA = Relative humidity field. 


ZZ1 = Height of cloud base. 


DZ = Vertical increment (meters) to be used. 


NCALL 


Number of times this subroutine called. Used as 


variable to bypass a READ statement. 


ieah 


li 
| 


coordinate for set of data coming in. 


JPT = J coordinate for set of data coming in. 


9. Calling arguments of GRPHCL and GRAPH are all parameters which 


have been computed or set in subroutine CLOUDS. 


10. Subroutine METMAP 


This subroutine was not written by the author, but is part 


of the Naval Postgraduate School Computer Facility Library and has 


been included for information. 


It is simply a shading routine that 


will print out a field of data with a 0.5 inch grid spacing and contour 


Ener Lie ld. 
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Y = Two-dimensional field to be contoured. 

N = Number of rows I in the array to be contoured. 

M = Number of columns J in the array to be contoured. 

T = Title for printout . Up to 96 columns of alpha-numeric 
information, 


BND = Bandwidth desired for contouring. 

AZ = Scaling constant; that is, each value will be multiplied 
by AZ before printing it out. Only the first three numbers 
to the right of the decimal point are on the output, so AZ 
can be used to control which numbers are to be printed out. 

BZ = Additive constant. This is useful when working with D-values. 
Usually, however, it is 0.0. If used, this constant will be 
added to each value of Y before printing it out. 

AMIN = Minimum value for subroutine to start contouring. 

IJT = 0 means subroutine will compute the minimum value and 

start contouring at that value. A value must still be 


specified for AMIN but it will not be used. 


ICON 1 if contouring desired. 


QO if no contouring desired. 

When working with this subroutine, it must be remembered that the 
grid point I=1, J = TIT, 1s the upper left corner of the grid, anda 
increases downward while J increases to the right. If an array is 
defined with the lower left corner as the point (1, 1), the follow- 
sat of FORTRAN statements will get the array in the proper order 


for the subroutine: 
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DO 5 I=1, 24 
DO 5 J=1, 24 
5 DUMMY (J,1)=Z(1,24-J+1) 


The field called DUMMY is then taken into METMAP as Y. 
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